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We compute the coefficient of restitution, e, starting from a microscopic model of elastic disks. 
qq ■ Collisions are found to be in general inelastic with a finite fraction of translational energy being 

| transfered to elastic vibrations. The coefficient of restitution depends on the relative velocity of the 

Q\ , colliding disks, such that collisions become increasingly more inelastic for larger relative velocities. 

t-H ' The predictions of the model are in agreement with the approach of Hertz in the quasistatic limit. 

The coefficient of restitution depends on the elastic constants of the material via Poisson's number. 

The elastic vibrations absorb kinetic energy more effectively for soft materials with low ratios of 
y-^- ' shear to bulk modulus. 

^ ! PACS numbers: 46.30.My, 62.30.+d, 83.70.Fn 
\0 ■ 
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, '■ I. INTRODUCTION 

> ■ 

The most important characteristic of granular media is the inelastic nature of the interparticle collisions. The 
I/"") \ removal of kinetic energy in a granular gas is responsible for nonequilibrium phenomena that are of theoretical and 
experimental interest. In computer simulations of granular matter this energy loss is usually treated in a very simplified 
way. In event driven simulations Q| a fixed coefficient of restitution is used, i.e. after each interparticle collision a 
certain fraction of the energy involved is lost. In experiment Q the coefficient of restitution is found to depend on 
velocity. When using molecular dynamics techniques || ad hoc phenomenological assumptions for the intergrain force 
, laws are introduced. (For a recent review, see Q.) 

Several microscopic mechanisms for the decay of kinetic energy during collisions have been discussed. Permanent 
plastic deformation of granular particles has been proposed as a possible mechanism for the removal of kinetic energy 
Viscoelastic behaviour was used to extend the theory of Hertz || to inelastic impact. One either invokes a 
phenomenological damping term in the equations of motion Q or uses a quasistatic approximation for low relative 
Ch , impact velocities Q. 

More recently temporary storage of energy in elastic modes fjjjLOll has been discussed for one-dimensional rods. 
During collisions energy is exchanged between translational motion and internal degrees of freedom, whereas in 
between collisions the energy stored in the elastic modes decays. If the collision rate is high enough, which is 
frequently observed in simulations as a precursor for inelastic collapse, elastic energy may be transformed back into 
kinetic energy. Thereby inelastic collapse is avoided and rich dynamics of temporary clusters is observed, including 
breakup and reformation of clusters as well as relaxation to a final state with all particles in contact jllj . 

The quantity controlling the fraction of kinetic energy, which is lost in a collision, is the coefficient of restitution 
e. We will restrict ourselves to head-on, normal collisions, where e is simply defined as the ratio of relative velocities 
after and before collision 

e := -V{/Vi . (1) 

For perfectly elastic one- dimensional rods phenomenological wave theory yields e = h/h, independent of the relative 
velocity of the colliding rods. Here l\ (I2) denotes the length of the shorter (longer) rod. If both rods are identical 
we have e = 1 always. Using the approach of Hertz, Rayleigh jl2| has estimated, that the fraction of energy stored 
in the fundamental mode for two slowly colliding spheres is about 0.02 • v\jc (where c is the velocity of sound in a 
one-dimensional rod of the same material). While being quite small for many realistic scenarios his value for the 
coefficient of restitution thus depends on velocity and does not vanish for identical spheres. 

In this paper we shall analyze collisions of two-dimensional clastic disks. Starting from three-dimensional clastic 
objects, the two-dimensional case can be realized in two ways, called plane stress and plane strain p5fl . Plane stress 
describes the situation of thin plates, where the stress components on the faces of the plates disappear. With the 
additional assumption of in-plane oscillations only, the equations of two-dimensional elasticity apply. In the case of 
plane strain the end sections of a prismatic body are confined between smooth rigid planes, so that displacements in 
the axial direction are prevented. If the forces do not vary along the length for symmetry reasons, it may be assumed 
that there is no axial displacement anywhere. This again reduces the problem to two dimensions. It is easy to show 
that the equations for plane strain are the same as those for plane stress, provided one makes the substitution 
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v^JL- (2) 

1 — V A 1 — V 

where E denotes Young's modulus and v Poisson's number. For plane stress Poisson's number v has the three 
dimensional value and hence is restricted to —1 < v < 1/2. With the above transformation this corresponds to 
— 1/2 < v < 1 for plane strain. In particular values of v close to 1 can be achieved for plane strain by materials with 
low shear modulus /i, like rubber. 

As in our previous work for one-dimensional elastic rods 10 , we model the transfer of translational energy 



to elastic vibrations of two-dimensional disks and discuss the collision between two identical disks or equivalcntly 
the collision of one disk with a rigid wall. The interaction is modelled by a hard core potential depending on the 
instantaneous separation between the two disks. In contrast to the one-dimensional case the equations of motion for 
the elastic vibrations have to be solved numerically. The main result of our paper is the coefficient of restitution as 
a function of initial relative velocity and elastic properties of the disks, as shown in Figure [l| Collisions are found to 
be elastic for vanishing relative velocity in agreement with Hertz' quasistatic approach || and become increasingly 
inelastic for increasing relative velocity. For soft matter, characterized by a small value of the shear modulus, the 
clastic modes provide a rather effective mechanism for the uptake of kinetic energy during collision, so that the 
coefficient of restitution decreases with decreasing shear modulus. In contrast to one dimension, where the collision 
of two rods of equal length is always perfectly elastic [pTof] , we find that the collision between two identical disks is in 
general inelastic. 

Our paper is organized as follows. In section p we discuss two-dimensional elasticity and calculate the eigenmodes 
of a disk with force free boundaries. In section HI we analyze the equations of motion of an elastic disk, colliding with 
a hard wall. The numerical solution of the equations of motion are discussed in section IV. In section ^ we apply 



Hertz' and Rayleigh's methods to two dimensions to study the limiting case of very small relative velocities, i.e. the 
quasistatic limit. We conclude with a summary of our results and an outlook. 



II. ELASTIC EXTENSIONAL MODES OF DISKS 



A. Elasticity of planar bodies 

We briefly review the theory of linear elasticity in two dimensions, as discussed e.g. by Love Landau jl4j or 
Timoshenko and Goodier [fil|. We consider small displacements u(x) and expand the free energy up to quadratic 
order in the strain field, which is defined as the symmetrized derivative of the displacement vector: 



1 / dui duj 



2 V dxj dxi 



(3) 



An isotropic elastic body has only two independent elastic constants, the shear modulus /i and the compressional 
modulus K. In terms of these the elastic free energy is given by 

d 2 x {n{ea - 5 t3 e u /2) 2 + Ke 2 u /2) . (4) 

The notation implies summation over indices which appear twice. The stress tensor is defined as the derivative of the 
clastic free energy with respect to the strain field and explicitly given by 

Vij = ^~f clah = 2^(e. y - - Sijeu/2) + K5 tj eu . (5) 
oeij 

It is sometimes convenient to introduce Young's modulus E and Poisson's number v instead of the moduli of torsion 
and compression. These quantities are simply related, in two dimensions we have E = (AfiK)/(K + /i) and v = 
(K — (j,)/ (K + /i) , so that alternatively to Eq. (^) we may use 

((1- v)eij + vSije u ). (6) 



lJ (l-^ 2 ) 

In the absence of body forces the equations of equilibrium read 

do 
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B. Normal modes 



In the following we shall discuss the eigenmodes of a circular disc with force free boundaries (see also [Eli). Our 
starting point are the linear equations of motion ]l4| ] 

V(V.u) + — Au=«L_^l_ 

(8) 

where p is the mass per unit area. The differential operator £ is hermitean for force free boundary conditions, as 
shown in the Appendix. To solve the coupled differential equations for u x and u y , we introduce the areal dilatation 
rj and rotation £ 

du x du v „ > du,, di^ 
Ox Oy Ox ay 

The equations of vibration (@) then simplify to 

Differentiating (|l^) with respect to x, ( pi] ) with respect to y and adding them yields 

_ 9 p(l — I/ 2 ) C? 2 77 , , 

^ =£L ^M- (12) 

Differentiating ( |l0| ) with respect to y, (11) with respect to x and subtracting them gives 

Hence, the equations of two-dimensional elasticity have been reduced to two scalar wave equations which are coupled 
through the boundary conditions. 

To discuss the eigenmodes of a circular disk, we use polar coordinates r, <f> with the origin at the center of the disk. 
The two-dimensional displacement vector is written as u = e^. The radial and tangential 

displacements u r and are given by 

u r — u x cos 4> + u y sin (f> , — —u x sin <f> + u y cos <f> . (14) 

Similarly strain, dilatation and rotation can be worked out in polar coordinates 

v = ^ + - + -^, 2C = r^ + ^_ir^ (15) 

or r r ocp ' ' 

du r u r 1 diii/, 

^rr = — ^ ; ^4><P — ' a~7~ > ^rch 

or r r dtp 

Finally, because the medium is isotropic the stress-strain relation in polar coordinates is as simple as it is for cartesian 
coordinates Eq. (^): 

Urr = Z J ^ rr Vt <l><P> i = i 2"( e, ^ + ^ err )' Cr ' , = 9 /-i T — T 6 ^' U^) 

The solutions to (|l^) and (|l3|) are of the form 

rj = A n K 2 J n (nr) cos n(f> cos , 2£ = B n K 12 J n {n'r) sinn0 cosui , (17) 
where J„ is Bessel's function of order n tlq], A„ and £?„ are constants and 
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(l-^ 2 )p 

E ' 



2(1 + x/)p 
E 



(18) 



There is a second set of solutions with "cos" and "sin" interchanged. We shall only consider situations which are 
symmetric in 0, so that u r (r, <fi) = u r (r, —cj>) and u^(r, cf>) — — uAr, —(p). Hence we restrict ourselves to the above set. 
The displacements corresponding to the above solutions are given by 



u r (r, cj>) 



A, 



dJ n (nr) 
1 dr 



nB n 
+ B. 



J n (/tV) 



r 

dJ„(re'r) 



dr 



cos n<f> cos cut , 
sinn<A cosuit . 



(19) 



This can be easily verified by substituting Eqs. (|l9| ) into Eqs. @, thereby one recovers the symmetric solution of 
Eq. (§). 

At the boundary of the disk, r — R, shear and compression have to vanish: a rr (R, (f) — and <r rr j ) {R 1 (f) = 0. For 
ri = and B = we have purely radial vibrations in which vanishes and u r is independent of <f>. The boundary 
condition then reads 



dJi(KR) 
dR 



R 



Ji(«i?) = 



(20) 



Areal dilatation and stress are maximal at the center of the disk, although the center is not displaced. The displacemnet 
vanishes along node lines, which are circles for n = 0. 

There is a second set of solutions with u r = 0, u$ independent of <f> and node lines which are cirles. These modes 
ar not excited in head-on collsions for which we impose the symmetry u,p{r, <fi) = —u^,{r 1 — 0). 

In the general case, i.e. for n > 1, there is no tangential displacement, = 0, for (/> = n an d no radial displacement, 

u r = for cf> = ^^n^ w ith k 6 {0 . . . 2ri — 1}. These values of represent n node lines for either or u r . The 
boundary conditions in general imply two equations 

l — i/ dJ n {nR) 



R dR 
+ nB n {l-v) 



1 



R? 

1 dJ n (K'R) 

R dR 



n 2 )J n (KR) 



1 . , ' 

•jjg J n (k R) 



(21) 



-2nA n 



1 dJ n (nR) 1 
R^R— & Jn{hR) 



Bn 



2 dJ n {n'R) 
R dR 



k' 2 - 



2r^ 
i? 2 



J„(K'i?) 



= 



We eliminate A n and i?„ and use 



to obtain the eigenvalue equation 



dJ n (nR) 



Kj n -i(K.R) 



Jn(nR) 

' R 



(22) 



(23) 



(1 — ^) 2 (1 — n 2 )nn J n _x{nR) J n _i(«;'i?) 
+ (1 - i/) (k 2 - (1 - i/)(l - n 2 )n) (Kj n ^ 1 (KR)J n {K , R) + k 1 J n -i(n' R)J n (nR)) 

+ (k 4 - 2(n 2 + n)(l - v)k 2 ) J n (KR)J n (n' R) = 



(24) 



For any fixed number n, there are infinitely many solutions u) n ,u numbered by I = 0, . . . , oo. These solutions have 
been determined numerically. Sometimes solutions for neighbouring I are very closely spaced and difficult to find 
numerically. It is then advantageous to simultaneously determine the zeros of the derivative of ( p4| ) , which are much 
more regularly spaced. In between a pair of zeros of the derivative, we then search for a zero of the function. 

The ratio A : B is determined by reinserting the values for ui n j into ((2^). We are left with one free constant for 
each eigenmode. This constant can be fixed by using normalized eigenmodes, according to 



dr r((t£''(r) cos(n0)) 2 + {u n /{r) sin(r^)) 2 ) = ttR 



(25) 
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Here we have introduced dimensionless quantities it"''(r) and u^' l (r) for the radial variation the displacement (Eq. [L9|) 

Ru?' l (r) := 



' dJ n (n nil r) JnW n ,i r ) 

A n ,l : h ™-Dn,Z 



i? U ^(r) :=- 



dr 

. JnWn,l r ) . D dJ„(K„,;r) 
™Ali h S n i 



dr 



Zeroes in it™' (r) and u.' (r) give rise to radial node lines of u r (r,(j)) and ^(r 



(26) 



A given value of / does not 



imply a given number of radial node lines of the radial or tangential displacement, because the boundary conditions 
mix the two displacement components. The displacements for a few eigenmodes are sketched in Figure ^[ and the 
frequencies, uj n ,i of the modes are plotted versus the number of radial nodes n in Figure ^[ For collisions, an important 
characteristic of a mode is the maximal radial and tangential displacement at the edge of the disk 



\R) 



(R) 



(27) 



The modes can roughly be divided into two classes with primarily radial or tangential displacement. In Figure y 
radial oscillations are characterized by C n j > S n j, tangential oscillations by C ni i < S n j. 

For n > 2 two distinct regimes can be distinguished. For small / the solutions are regularly spaced and the 
arrangement barely changes, when going from n — > n + 1 . The difference in magnitude between tangential and radial 
displacement is small (e.g. C„,i — S n j), and radial as well as tangential modes tend to cluster, the more so the larger 
n. For I 3> n the solutions, w„j for the radial and tangential modes appear to be almost independent of each other, 
as they are in the case n — 0. This results in an irregular spacing of the frequencies of the modes. The distinction 
between primarily radial and tangential displacement is pronounced. Ratios C n j : S n ,i typically are larger than 1 : 10 
(or 10 : 1) when I 3> n, exceeding 1 : 100 for very large I. When two solutions have very similar frequencies the ratios 
are smaller at the surface, but this does not change the overall character of the mode. The number of node lines 
increases uniformly with 1 in both cases. As a general statement we conclude, that the mixing of the displacements 
introduced by the boundary conditions becomes more important with increasing n and decreasing I. 

For n = 1 dilatation and rotation, stress and shear disappear linearily at the origin, but there is a finite displacement 
of the center (see Figure ^). For n > 2 displacement of the center disappears like r n_1 , dilatation and rotation like 
r™, stress and shear go like r™ -2 (and are maximal at the origin for n = 2). Thus for increasing n the deformations 
are increasingly concentrated towards the edge of the disk. 

The mode with the lowest frequency is the quadrupolar mode (n = 2, 1 = 0), where contraction in one direction is 
met by expansion in the orthogonal direction (see Figure ||) . This mode will turn out to be the most important mode 
for the storage of elastic energy, as suggested by Rayleigh for three-dimensional spheres fll2j. 



C. Energy of vibrations 

We have a complete orthonormal system of eigenfunctions such that any (symmetric) state of the disk can be 
expanded in this set: 

u r (r,<f>)\ _ ^ ( u?> l (r) cos (ncf)) \ 

U0 (r» ) ~ IJ 4 ^ \ u n /(r) S m(ncj>) ) ' ^ j 

To express the elastic energy E e \ as in terms of the expansion coefficients, Q n ,i, we first use partial integration together 
with force free boundary conditions to rewrite the elastic energy as 

E f 

E c i as = -— / d 2 xVu fe (x)(£u) fc (x) (29) 

2(1- v 2 ) J s 

Next, we use the equation for the eigenfunctions of C 

(£u"' i ) fc (x)=-4 ii <' i (x) (30) 

to rewrite the elastic energy as 

-Solas = -Q-y^n.aQnJ- ( 31 ) 

n.l 
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The kinetic energy of a vibrating disk is given by 

2tt r 

Skin - | / ck/> y dr r(u*(r, 0) +«J(r,0)) (32) 



so that the total energy of a vibrating disk not interacting with external forces is given by 

H ° = E + ' (34) 

where we have introduced canonical momenta P n j — mQ n j. In the elastic approximation a vibrating disk behaves like 
a set of independent harmonic oscillators with canonical loci Q n j and canonical momenta P n ,i, obeying Hamilton's 
equation of motion 

dQ„, ; 0H° P n>/ AP ntl dH° 

at aP n ,i m at oQ n ,i 



III. EQUATIONS OF MOTION IN THE PRESENCE OF A WALL 

The head-on collision of two identical elastic disks is equivalent to one elastic disk hitting a hard wall, which itself 
is not deformed elastically. We choose a coordinate system such that the wall is located at x — (see Figure [|) and 
the disk is approaching from the left. The interaction between disk and wall is modelled by a potential V(x) — Voe ax , 
where x = x(4>) is the distance between the edge of the deformed disk and the y-axis. The choice of Vo is arbitrary, 
because it only affects the position of the wall, as can be seen from the substitution Voe ax = y^ e a ( x +{^ k )/ a ) with 
Vb = /cVq. A hard core interaction is recovered in the limit a — > oo. We assume that the disk is moving along the 
a;-axis in positive direction. Its center of mass position is denoted by xq (t) . 

We expand an arbitrary, symmetric state of the disk in normal modes 



u r (r,(f>,t) \ _ sr^ n , K' l (r)cos{n(/))\ 
ut(r,<j>,t) J ~ ^ n ' lW \u n /(r)sm(ncl>) J 



(36) 



with time dependent expansion coefficients Q n ,i(t)- Since the wall is only pushing against the edge of the disk we 
need to know the location of the boundary. The radial and tangential distortion at the edge are given by 

u r {R,(j>,t) = ) j Q n ,i(t)C nt i c °s n<t>, u$(R,(j>,t) = )^ Q n ,i(t)S n j sinnift . (37) 

n,l n,l 

The location ( ,^'5 | of an element at ( ^o(*)+ -Rcos \ ^^^^ d e f ornia ti on is then given by 
V 2/OM) / V Rsmcf) J 

x(4>,t)\ _ ( x (t) + (R + u r (R, <f>)) cos<f> — u<f,(R, 4>)sm.(f) \ , . 

y{(f>,t) ) ~ \ (R + u r (R,ct))) sin cjj + u^R^) cos <j) )' 

Its rr-component enters into the wall potential and is expressed in terms of normal modes as follows 

x(<f),t) := xq (t) + R cos tf> + V^ Q n i (t) (C n> i cos n(j) cos <f> — S n> i sin ruf) sin 4>) . (39) 

n . / 

The total energy of an elastic disk interacting with a fixed wall is given by 

tt/2 

H --=^Po + f2{^ P li+ m <iQli)+ V o J dct>e™W. (40) 

n ' 1 -tt/2 
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The dynamic evolution of the expansion coefficients Q n ,i(t) follows from Hamilton's equations of motion 
dQ nJ = dH_ 

at dP n ,i ' 

dP, hl _ dH 



dt dQ n ,i 

tt/2 

— muj nt iQ nt i + aVo J d(f){C n ,i cos 770 cos</> — S n ,i sin77</> sin</j) & ax ^-^ . (41) 

-7T/2 

Instead of using real valued functions Q n ,i{t) and P n ^(t), we find it convenient to introduce a complex function 
q n ,i(t) = q^i(t) + iqi t i(t) such that 

Qn,l{t) =Rc(g„,,(t)e iu "'' t ) , 

P„,,(i) = mLu n jIm(q n j{t)e iu ^ . (42) 
Hamilton's equation of motion then read 

4n,l C0S (Un,lt) ~ qi t l Sin (u n jt) = , 

tt/2 

q^ l sin (u> n it) + ql , cos (u) n it) = - - / d<f>(C n ,i cos n<J} cos </> - 5„,; sin 770 sin </>) e ax ^-^ . 

-7T/2 

These two equations can be combined into a single equation for the complex function q rh i{t) 

tt/2 

d q n j = -l^-e^.it [ d0(C n , z cosn</) cos0- 5„jsinn</> sin0) ■ e ax ^'^ . (43) 



dt ' niuj n j 

-tt/2 

The time evolution of the center of mass velocity also follows from Hamilton's equation of motion 

d« = j_d|o = _i_ e H = _?Vo f d ^ )t) . (44) 

dt 777 dt 777 oxo m J 

-7T/2 

completing our description of a two-dimensional elastic disk hitting a wall. The only values needed for every mode 
(n,l) are the frequency ui n ,i and the maximal displacements C n ,i and S n ,i- The actual numerical integration of the 
equation of evolution (^3|) requires the consideration of many different modes to find a good estimate of the coefficient 
of restitution. 



IV. NUMERICAL RESULTS 

The actual integration of (^3|) is numerically demanding, when high accuracy is required. Time steps have to be set 
quite low to capture the oscillations of the fastest modes. Numerical accuracy is limited by the number of modes that 
can be used (usually less than 1000), so that a 4th order Runge-Kutta method suffices for the numerical integration 
of the equations of motion. Starting from x(4>,t) and q n ,i{t), we update xq according to Eq. ( f44| ) and all n n j for 
a given set (77 < 77 max ,Z < Z max ) according to Eq. ((43|). Then a new x (<M + At) is calculated from Eq. rt3|). For 
a normal collision of a disk with no excited modes the integration over the edge of the circle can be limited to the 
upper quadrant because of symmetry. We typically discretize the upper quadrant with 256 values for </> and retrieve 
the values of cos mf> and sinruj) from a look-up-table to effectively speed up the computation. Initial conditions are 
Qn.i = f° r all n, I, and a value of Xq sufficiently negative, such that there is no interaction with the wall. 

In our simulations and plots we measure velocity in units of c = W E / p, distances in units of R and forces in units 
of mc 2 /R = ttRE. The numerical simulations have been performed for a range of initial velocities 0.005c < v; < 0.3c 
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rank 


n 


1 


average energy 


energy after collision 


1 


2 





0.1594843 


0.0513636 


2 


3 





0.0480930 


0.0234695 


3 








0.0231257 


0.0089346 


4 


4 





0.0220028 


0.0038760 


5 





1 


0.0128242 


0.0003144 


6 


5 





0.0114396 


0.0027244 
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1 


1 


0.0110218 


0.0002610 


8 


1 





0.0066413 


0.0002548 


9 


6 





0.0059540 


0.0021365 


10 





2 


0.0048560 


0.0000279 



TABLE L The ten most important modes ranked according to the average fraction of the total energy during collision for 
a = 500-R, Vi = 0.1c using 1600 modes. 

and a a value of v = 0.33 for Poisson's ratio. A few additional runs were made with other values of v to check the 
dependence of our results on Poisson's ratio. In Figure [| the collsion of an elastic disk with a rigid wall is shown, as 
calculated numerically according to the above procedure. 

Computation of the coefficient of restitution requires two different extrapolations. Modelling disks as elastic continua 
requires the limit of infinitely many modes. At the same time the artificial exponential potential, introduced to make 
the collision accessible to numerical calculations, should be made infinitely hard. 

A harder potential of the wall increases the energy, which remains in the elastic modes for two reasons: 1) The force 
which excites the elastic modes is stronger. 2) During the collision less energy is stored in potential energy, which 
after the collision is returned to kinetic energy and hence cannot be transferrred to the internal degrees of freedom. 
Raising the number of modes has the opposite effect. The coefficient of restitution e is closer to one for a larger 
number of modes, because the larger number of modes allows a better adjustment of the shape of the disk, so that 
collisions become effectively "softer" . In both cases extrapolation from finite values is possible, it is more difficult for 
the number of modes than for the potential parameter a. 

To select the right set of modes, we use the following procedure. For every initial velocity v\ the modes were ranked 
according to their average energy content in a preliminary run. As expected, the average energy content of a mode 
decreases with increasing frequency w n> i. It turns out that the purely radial modes (see Eq. (|2^)) are particularly 
important. Typically 75 % of the most significant modes have n = 0. Another essential class of modes has 1 = 0. 
The importance of these two classes can be inferred from Figure ||, where we show the distribution of energy over 
modes (n,l) at closest appraoch. In Table 1 we list the energy averaged over the entire collision as well as the energy 
after collision for the ten most important modes. It can be seen that the average energy content of a mode during 
the collision can be very different from the energy remaining in the oscillation after the collision has been completed. 

In Figure [?| the energy which is absorbed in the elastic modes is plotted versus potential parameter a. It scales 
nicely with l/a, because the energy that is stored in the wall potential at closest approach is also oc l/a. Hence 
the extrapolation a — > oo is straightforward, except for small initial velocity, when the coefficient of restitution for a 
given potential parameter a approaches one. To identify the inelastic signature of the collisions one has to choose ever 
harder potentials. This makes it difficult to extrapolate to the limit Vy — > numerically. Fortunately the quasistatic 
limit can be treated analytically, as discussed in the next section. 

The extrapolation of the energy stored in vibrations to infinitely many modes is shown in Figure |^. The limit 
N — > oo is approached like 1/yN. The more violent collisions need a higher number of modes for the 1/y/N regime 
to appear. 

The two extrapolation procedures have been performed for a range of initial velocities and two values of Poisson's 
number. The resulting coefficient of restitution e is plotted versus the initial velocity in Figure [l]. It is seen to approach 
one in the quasielastic limit and to decrease with increasing relative velocity. Collisions are found to be more inelastic 
for higher values of v. 

V. QUASISTATIC APPROACH 

As we just have seen, the limiting case v-Jc <C 1 is difficult to treat numerically. However Hertz' law of contact 
can be extended to two dimensions. Solving the problem in two dimensions is actually more complicated than the 
three-dimensional case, because there are no local solutions. 



8 



The quasistatic acceleration of a disk by a hard wall is equivalent to a disk in equilibrium in a gravitational field 
and supported by a hard wall. This problem has been solved for a point contact by H. Michell 0|. To compute the 
compression caused by the gravitational field we have to generalize his solution to an extended contact between the 
disk and the supporting wall. 

In the contact pressure is calculated for a two-dimensional contact of cylindrical bodies. This solution can 
be taken over to our problem of a two-dimensional disk with one-dimensional loading. Inside the loaded region 
— a < y < a (see Figure ^) a normal stress p(y) given by 

P(y) = ^(a 2 -y 2 ) 1/2 (45) 

acts on the surface. Here P is the total load, R is the radius of the disk and a 2 = APR/ (irE). Integrating over the 
loaded region yields the stresses in the elastic disk due to the normal pressure p(y): 



(Tyy(X,y) 

cr xx (x,y) 



2x f p(s)(y — s) 2 ds 



:,,r - ,: - I [{y-s) 2 + x 2 } 2 

— a 

a 

2x 3 f p(s)ds 



n J [{y — s ) 2 + x 2 } 2 



2x 2 f p(s)(y — s)ds 



t J [(y- 



s) 2 



(46) 



Along the axis of symmetry integration is possible and one finds 



(j yy (x, 0) = ( ° , + ^ 1/2 - 2a;) , i 47 ! 



ira 
2P 



2)1/2 
,2 , _2\-l/2 



o- M (x,0) = (a'+x^) x ^ , <T«„(aj, 0) = 0. 

In a three-dimensional elastic medium, the displacement decreases like 1/r for large distances r from the loaded 
region. In a two-dimensional elastic medium with one-dimensional loading the displacement varies as ln(r), so that 
the displacement can only be defined relative to an arbitrarily chosen datum (which for three-dimensional systems 
is usually taken at infinity). This implies that the total deformation 8 cannot be computed from the local stress 
distribution p{y) . To find the total compression for a two-dimensional system one has to consider the stress distribution 
in the bulk of the two-dimensional elastic body as well as its shape and size. 

To this end we need to generalize the solution of Michell |l7j , who solves the problem of a heavy disk supported by 
a point force. The stress generated by the gravitational force is most easily evaluated in a coordinate system whose 
origin is located in the center of the sphere. The equations of equilibrium for a uniform downward force of magnitude 
pg read 

da xx da xy da yy da xy 2 

-^ + ^=P9, ^ + ^T=°> VK,+^) = (48) 

and are solved by 

1 1 1 

°xx = -^pgx , a yy = --^pgx , a xy = -pgy . (49) 

The stress due to the point contact is purely radial in a coordinate system whose origin lies in the point of contact 

o cos <h 

a rr = -2pgR 2 —— . (50) 

To obtain the total stress distribution we transform the stress due to gravitation ( |49| ) to the coordinate system whose 
origin lies in the point of contact 



^pg(x - R) , a V y = -^pg(x - R) , a xy = ^pgy. (51) 
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and superimpose the two contributions to the stress. The resulting solution does not satisfy the boundary conditions 
of stress free edges, which can be guaranteed by adding a uniform biaxial tension a xx = a yy = pgR. Hence the total 
stress is given by 



pg(x/2- 



2R A x 



2^3 



-P9 



(x 2 + y 2 ) 2 . 

( (l -2 fl )/2 + -^). (52) 



The radial stress distribution ( j50| ) from the point contact (which creates a logarithmic divergence in the displace- 
ment) will now be replaced by the realistic stress distribution (|47|) for an extended contact. The total stress along the 
axis of symmetry y — is then given by 

<T XX (X,0) = --( (fl2+ 3.2)1/2 " 2^) ' 

P ( 2(a 2 + 2x 2 ) 4x x- 2R\ . . 

^(^.0) = -- I 9.9 , -9M/9 + on2 ' ( 53 ) 



it \a 2 {a 2 +X 2 ) 1 / 2 a 2 2R 2 J 



Here the total load is simply given by P — mg = irR 2 pg. 
The strain can be calculated from the formula 

t-XX = ^{VXX - VOyy) (54) 

hi 

and the total compression S of the disk is found by integrating e xx from x = to x = 2R. For a <g; R we find 

2R 

e xx dx = -?-(2 ln(4i?/a) - 1 - v) (55) 


P /, AR-kE. \ , , 

= ^(m{— )-!-»)■ (56) 

The compression S can also be obtained from the related problem of a disk compressed between two walls, which 
has been solved in [^9|. In this case the total deformation is twice the deformation as calculated above for small S. 
In Figure |l^ we show P(S) as computed by three different methods. The analytical results are compared to the 



dynamic calculations of Section IV. If the disk is held fixed at a given distance from the wall, it will assume a 
deformation which minimizes its total energy. This deformation and the force actingon the disk can be computed 
using the expansion of the elastic deformation by the modes calculated in section y. The total energy given by 
Eq. (EOT) with P n _i = is iteratively minimized until convergence is achieved. The numerical result of this calculation 



is shown as a dashed-dotted line in Figure 10, The results from the quasistatic calculation and from the minimization 
procedure agree well for small S, whereas for larger deformations one observes deviations from Hertz law of contact. 
The forces during the dynamical collision are considerably higher than in the quasistatic limit, because the compression 
is localized. The results of the dynamic calculation are shown for two initial velocities to demonstrate the approach 
towards the quasistatic case for the smaller initial velocity. 
In the limit of small P, Eq. (B6|) reads 



and may be inverted to yield 



P , ARttE, 

6 = ^(—), (57) 



SttE ,<nnlnlM, , x 

ln4R/5 v (ln£) 2 ' 



To lowest nontrivial order the potential energy is given by 

^« ■ < 59 » 

The kinetic energy before collision is given by ^vfirR 2 p. In the limit of low Vi, where the quasistatic approximation 
is expected to hold, conservation of energy can be used to obtain the maximum value of 5 
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y- / 4c R / 4c 

{ max ~ —R\ In — and hence r cx — Win — (60) 
c V Vi c \ Vi 

for the contact time r. In the limit t> i — > the contact time diverges logarithmically, i.e. much slower than for three- 
dimensional spheres, where Hertz' theory predicts r oc v i 1 ^ 5 . In Figure [ll] contact time and maximal compression 
are plotted with the force calculated from Eq. (|5^) as a function of the initial velocities. The very slow divergence 
of the contact time with m — > can be observed. The results for the contact time as obtained from the quasistatic 
calculation agree well with the numerical simulations of the full dynamic problem in the range of velocities were both 
methods apply. 

Our aim is an approximate expression for the coefficient of restitution e in the limit of small velocities Wj. This 
can be achieved by using the quasistatic force law in the equations of motion. Our approach is a generalization of 
Rayleigh's [ fL2| , who derived the energy stored in the fundamental mode in the case of spheres. 

We replace the wall potential in Q4Q ) by a potential V, that acts only at the point <j) = of the boundary. This 
amounts to 

1 °° 1 

H := 2^Po + E ( 2^ P li + m <iQh) + V ^ = °' *)) > (61) 

n,l 

where x{4> = 0, t) = x Q (t) + R . + ^ Q n ,i(t)C n ,i . (62) 

n , I 

From Hamilton's equations of motion we derive 

77 = TT77 = mUnlQnl + G„.i— . (63) 

dt oQnj dx 
Next we approximate dV/dx by the force P as calculated from (|56|). Using the compact notation ( p2| ) we can write 

-^ n ,it_^n£_pU) , (64) 

P(t) , (65) 

(66) 

These equations can be integrated numerically with initial conditions q n j = 0,5 = 0,v = V{. The energy which 
remains in the modes after collision is given by 

H cx = — ^ ^ n ,Mn,i{ T )\ 2 ■ ( 67 ) 

In Figure |l2| we show how much energy is stored in the modes for low V{. The quasistatic approximation agrees in 
magnitude with the result from the dynamical calculation. Most of the energy goes into the quadrupolar mode. Since 
the quasistatic approach is only correct to first order, the decreasing vibrational energy for Vi > 10 _2 c is presumably 
unphysical. 



dt 


— le 


d 


1 


~r v — 




dt 


m 


±5 = 






—v 


dt 



VI. CONCLUSIONS 



It was our aim to derive the coefficient of restitution as a function of velocity, starting from a microscopic model of 
simple objects with internal degrees of freedom, which can absorb part of the kinetic energy of translation. We have 
discussed in detail the head-on collision of two clastic disks, initially nonvibrating. This problem is equivalent to the 
collision of an elastic disk with a rigid wall, representing the plane of reflection symmetry of the two colliding disks. 

The dynamics of a collision has been formulated with help of Hamilton's equations of motion for the normal 
coordinates of the elastic disk, interacting with a repulsive wall potential. The resulting dynamic equations were 
solved numerically for a finite number of modes and then extrapolated to the case of infinitely many modes. 

The main results are the following. In contrast to the quasistatic theory of Hertz, the dynamic collision of two 
identical elastic spheres is inelastic in the sense that the relative velocity is decreased, corresponding to a coefficient 
of restitution smaller than one. The amount of translational energy which is converted to vibrational energy depends 



11 



on relative velocity and Poisson's number v. The conversion is more effective for materials with low shear modulus, 
corresponding to a value of v close to one. For high initial velocities the coefficent of restitution can be rather small, 
e.g. e ~ 0.65 for v — 0.9 and i»i = 0.3c. In the quasistatic limit, V\ — > 0, collisions become more and more elastic. We 
have generalized the quasistatic approach of Hertz and Rayleigh to two dimensions and shown that the contact time 
r diverges logarithmically as i?j — » and the maximum deformation vanishes as 5 max ~ vit. The dynamic approach 
agrees with the quasistatic theory in the regime, where both theories apply. 

We expect the collisions to be more strongly inelastic, if the two spheres have different size. For one-dimensional 
elastic rods we could show that the coefficient of restitution is strictly one, if the two rods have equal length, no 
matter whether vibrations are excited before collsion or not. In the latter case, when no vibrations are excited, the 
coefficient of restitution is given by the ratio of lenghts of the two rods. We are presently extending our calculations 
to disks of different radii. Future perspectives include collisions of elastic spheres as well as more elaborate models of 
contact, e.g. including roughness. 
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APPENDIX 



In this appendix we show that the differential operator C as defined in (g) is hermitean with force free boundary 
conditions, which read 



<r rr (R, <j>\ 
a r ,p(R, 4>] 

We have for two arbitray displacement fields u and v 



1 1 + v , . 1 — v . 
v • ( — -— V(V • u) + — — Au 



8u r 


(U r 

Yv — + 
V r 


1 <9it0 


dr 


r dcj) . 




+ 1 9 




—<t> 

du r 


r dutj, 


r 



o, 



r=R 



r=R 



(68) 



da; 2 u- f — t^V(V ■ v) 



+ boundary terms 



1 - v 



-Av 



(69) 



C is hermitian, if the boundary terms vanish with u and v satisfying Eq. (|6q). Using the normal vector f along the 
boundary and ds = r dxj), the boundary terms read 

l + v 



J ds((V • u)(f • v) - (V • v)(r • V ; 

J d.s(v (r ■ V) • u) - u (f • V) ■ v; 



dS 
1 - v 



(70) 



l + v I" ( . du r u r 1 du^ , . dv r v r 1 dv<p 
2 J \ r dr r r dcj) r dr r r dcj) 

OS 



1 - v 



, du r dud 

dS ( V r — h Vrf, — U r — h U„ 



i(v r 



OS 



dr 

u r 1 du 



dr 



dv r 
dr 



dv4> 
dr 



(71) 



- [ ds(u f— +v(— + l^i. 
J \ r \ dr \ r r dcj> 

dS 



1-vf dvs 1 dv^ \ 

Mr— + U 4>--^T ) 



d<p 



1 - v / dud 



1 dUr, 



V6-- 

r 



(72) 



In (jT^) we rearranged the terms for our purposes. The following relationships are derived using partial integration 
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271 



2tt 



d(f> v^-^-j- = — I d(f> u 



o 

2tt 



d<p v^—— 



o 

271 




d(j) u^—— 



(73) 



Both times the boundary terms vanish because u and v are single valued functions of (j>. The same relationships 
hold with u and v interchanged. We use these and subtract a term u^v^/r from both integrals to finally write the 
boundary terms as 



, r (du r (u r ldus\\ l — v (Ou r 



1 — v / du r 1 du^ 



d(j) r d(j) 



} 



i)S 



/, f / dv r ( v r 1 dv$ \ \ 1 - v ( dv r 1 dv^ V<h \ i I 

y ds rK 97 + <7 + r ^) ) + — M n ^ + ; 90 " v) I = 

95 



(74) 



Obviously, the lefthand side of (|74|) vanishes if the boundary conditions (|68|) are fulfilled by u and v. Finally we 
observe that C is also hermitian with respect to a fixed boundary u r = = and v r — — 0. 
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Captions 

Figure [l| Coefficient of restitution e as a functio n of initial velocity for v = 0.33 and v = 0.9. Here u; denotes the 
initial relative velocity, expressed in units of c = J E / p = 1 (see below) . 

Figure ||: Plot of the displacement of three simple eigenmodes. Radii and circles of the originally undeformed disk 
assume the distorted shapes sketched here. 

Figure ||: Dimensionless eigenfrequencies ui n .i ■ R/c of an elastic disk versus n (v — 0.33). We distinguish radial and 
tangential oscillations, according to the dominant displacement at the edge of the disk. 

Figure H: Illustration of the model of an elastic disk hitting a rigid wall. 

Figure ]E[ Time sequence of a disk with v\ — 0.2c hitting a rigid wall located at the origin. The deformations 
remainingafter the the collision can be clearly identified. 

Figure |6j: Logarithm of energy (arbitrary units) in modes with low [n, I) at closest approach. The spectrum is 
dominated by the quadrupolar mode (n = 2, 1 = 0). 

Figure 0: Fraction of the energy stored in the clastic modes 1 — e 2 after collision with a wall potential V = Voe ax , 
extrapolated to a — > oo. The different curves (from bottom to top) correspond to different values of initial velocity: 
v-Jc = 0.005, 0.01, 0.02, 0.04, 0.06, 0.30. For U; < 0.1c only the three most significant data points were used to 
determine the regression lines. N = 940 modes were taken into account. 

Figure |[ Fraction of energy stored in vibrational modes extrapolated to N — > oo. The energy scales with l/*/N. 
For the regression lines N = 300, 400, 600 and 940 were used. 

Figure pj: Visualization of the quasistatic approach. 

Figure ]10| : Dimensionless force versus compression as calculated with different approaches. For the dynamic cal- 
culation we have used a/R — 500 and two values of initial relative velocity, v\jc = 0.1 and v-Jc = 0.04. The arrows 
indicate the direction of time. 1300 modes were used for the relaxation towards the minimal energy for a given 
distance. 

Figure [H]: Dimensionless contact time r ■ R/c and rescaled maximal compression (<5 max / R)(c/vi) as calculated from 
the quasistatic approximation. The result for r compares well with data points (O) from the dynamical calculation. 

Figure [l2|: Fraction of energy stored in the lowest 1800 modes versus velocity using the Rayleigh approach. For 
these very low velocites the energy uptake is dominated by the quadrupolar mode (n=2,l=0). For comparison we 
show data points (O) from the full dynamic calculation. 
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FIG. 1. 
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FIG. 2. 
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FIG. 10. 
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